1 Introducción

En este trabajo se analiza la estructura multicapa de la red de transporte de la ciudad de Madrid a partir de dos conjuntos de datos reales proporcionados:

Un fichero con las paradas y estaciones de distintos modos de transporte.
Un fichero con información de viajes mensuales por distrito.

Los objetivos concretos son:

Construir una red multicapa donde cada capa representa un modo de transporte (bus, metro, bicimad, parking).
Estudiar la centralidad de las zonas de intercambio en cada capa.
Medir la versatilidad de las zonas (cómo cambia su importancia entre modos).
Realizar una visualización multicapa en la que se vean explícitamente las distintas capas conectadas.
Explorar de forma sencilla la vulnerabilidad del sistema ante fallos en una de las capas.

Toda la parte empírica se basa exclusivamente en los ficheros:

transporte_madrid_consolidado.csv
viajes_madrid.csv

No se utilizan ejemplos externos (como StarWars) para que el análisis esté completamente centrado en los datos de Madrid.

2 Carga y exploración de datos

library(tidyverse)
library(janitor)
library(sf)
library(dbscan)
library(FNN)
library(igraph)
library(scales)
library(knitr)

2.1 Lectura de los ficheros

Se asume que los ficheros CSV están en el mismo directorio que este documento RMarkdown. Si no es así, ajustar las rutas.

ruta_transporte <- "transporte_madrid_consolidado.csv"
ruta_viajes     <- "viajes_madrid.csv"

transporte_raw <- readr::read_delim(
  ruta_transporte,
  delim = ";",
  col_types = cols()
)

viajes_raw <- readr::read_delim(
  ruta_viajes,
  delim = ";",
  col_types = cols()
)

glimpse(transporte_raw)

## Rows: 14,188
## Columns: 6
## $ stop_name         <chr> "                                                   …
## $ stop_lat          <dbl> 40.37835, 40.38341, 40.36813, 40.37773, 40.40753, 40…
## $ stop_lon          <dbl> -3.75458, -3.62481, -3.62283, -3.67311, -3.59176, -3…
## $ transport_mode    <chr> "bicimad                                            …
## $ stop_id           <chr> "382 -Parque Eugenia de Montijo, 2                  …
## $ distrito_asignado <dbl> 2807911, 2807918, 2807918, 2807913, 2807919, 2807911…

glimpse(viajes_raw)

## Rows: 8,064
## Columns: 11
## $ fecha             <chr> "2024-01", "2024-01", "2024-01", "2024-01", "2024-01…
## $ name              <chr> "Centro                                            "…
## $ edad              <chr> "0-25", "0-25", "0-25", "0-25", "0-25", "0-25", "0-2…
## $ sexo              <chr> "hombre", "hombre", "hombre", "hombre", "mujer", "mu…
## $ numero_viajes     <chr> "0", "1", "2", "2+", "0", "1", "2", "2+", "0", "1", …
## $ personas          <dbl> 90541756, 18809207, 73248991, 187551534, 75357987, 1…
## $ Codigo            <dbl> 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, 28, …
## $ Provincia         <chr> "Madrid                ", "Madrid                ", …
## $ distrito          <chr> "2807901 ", "2807901 ", "2807901 ", "2807901 ", "280…
## $ zona_pernoctacion <dbl> 2807901, 2807901, 2807901, 2807901, 2807901, 2807901…
## $ poblacion         <dbl> 140644, 140644, 140644, 140644, 140644, 140644, 1406…

2.2 Descripción básica de `transporte_madrid_consolidado`

transporte <- transporte_raw %>% 
  clean_names() %>% 
  mutate(
    transport_mode = str_trim(transport_mode),
    transport_mode = str_extract(transport_mode, "\\w+")
  )

table(transporte$transport_mode)

## 
## bicimad     bus   metro parking 
##     626   12430    1050      82

summary(select(transporte, stop_lat, stop_lon, distrito_asignado))

##     stop_lat        stop_lon      distrito_asignado
##  Min.   :40.28   Min.   :-3.875   Min.   :2807901  
##  1st Qu.:40.40   1st Qu.:-3.711   1st Qu.:2807905  
##  Median :40.43   Median :-3.689   Median :2807909  
##  Mean   :40.43   Mean   :-3.683   Mean   :2807910  
##  3rd Qu.:40.46   3rd Qu.:-3.656   3rd Qu.:2807915  
##  Max.   :40.56   Max.   :-3.447   Max.   :2807921

Variables clave:

stop_name: nombre de la parada/estación.
stop_lat, stop_lon: coordenadas geográficas.
transport_mode: modo de transporte (bus, metro, bicimad, parking).
stop_id: identificador de parada.
distrito_asignado: código de distrito (INE).

Esto define los nodos base del sistema de transporte.

2.3 Descripción básica de `viajes_madrid`

library(readr)
library(dplyr)
library(janitor)
library(stringr)

# Ajusta la ruta a tu fichero si hace falta
viajes <- read_delim("viajes_madrid.csv",
                     delim = ";",
                     locale = locale(decimal_mark = ","),
                     show_col_types = FALSE) %>% 
  clean_names()

# Comprobamos que personas es numérica y razonable
summary(viajes$personas)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    4515   37486  110989  140469  186820  737342

Variables clave:

fecha: mes de referencia.
name: nombre del distrito.
edad, sexo: características sociodemográficas.
numero_viajes: categoría de nº de viajes (0, 1, 2, 2+).
personas: número de personas en esa celda.
distrito: código de distrito.
poblacion: población total del distrito.

2.3.1 Indicador de intensidad de viajes por distrito

Construimos un indicador sencillo de viajes medios por persona en cada distrito. Primero pasamos la categoría numero_viajes a un valor numérico aproximado y luego calculamos una media ponderada por el número de personas.

2.3.2 Indicador de viajes medios diarios por persona y distrito

La variable personas se encontraba codificada con formato decimal europeo, por lo que fue necesario convertirla explícitamente a formato numérico. Asimismo, la variable numero_viajes es categórica (0, 1, 2, 2+), por lo que se transformó en una aproximación numérica asignando el valor 3 a la categoría “2+”. El indicador final de viajes medios por persona se calculó como una media ponderada por el número de personas en cada categoría, obteniéndose valores en torno a 2 viajes diarios por persona, coherentes con patrones conocidos de movilidad urbana.

library(dplyr)
library(knitr)

viajes_distrito <- viajes %>% 
  mutate(
    # Pasamos las categorías de nº de viajes a un valor numérico aproximado
    viajes_cat = dplyr::case_when(
      numero_viajes == "0"  ~ 0,
      numero_viajes == "1"  ~ 1,
      numero_viajes == "2"  ~ 2,
      numero_viajes == "2+" ~ 3,
      TRUE ~ NA_real_
    )
  ) %>% 
  group_by(distrito) %>% 
  summarise(
    # Total de personas (peso) en el distrito
    personas_totales = sum(personas, na.rm = TRUE),
    # Total de viajes ponderados
    viajes_totales   = sum(viajes_cat * personas, na.rm = TRUE),
    # Viajes medios por persona
    viajes_medios_por_persona = viajes_totales / personas_totales,
    .groups = "drop"
  )

kable(
  head(viajes_distrito),
  digits = 3,
  caption = "Viajes medios diarios por persona y distrito"
)

Viajes medios diarios por persona y distrito
distrito	personas_totales	viajes_totales	viajes_medios_por_persona
2807901	51888948	110673624	2.133
2807902	51542198	108276834	2.101
2807903	39251902	78164164	1.991
2807904	49039609	98112207	2.001
2807905	48349024	95741527	1.980
2807906	53679153	110757611	2.063

summary(viajes_distrito$viajes_medios_por_persona)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   1.956   2.004   2.038   2.041   2.063   2.135

El indicador de viajes medios diarios por persona presenta valores muy homogéneos entre distritos, situándose en un rango estrecho comprendido aproximadamente entre 1.98 y 2.13 viajes diarios por persona. Esta baja dispersión sugiere que, a escala distrital, los patrones de movilidad cotidiana en Madrid son relativamente uniformes, independientemente del volumen total de población o del número absoluto de desplazamientos. En particular, los distritos con mayor intensidad de viajes no necesariamente coinciden con los más poblados, lo que indica que la demanda de movilidad responde más a la estructura funcional del territorio y a la conectividad del sistema de transporte que al tamaño demográfico. Este resultado refuerza la idea de que la versatilidad multimodal y el papel estructural de ciertas zonas en la red pueden desempeñar un papel más relevante que la mera intensidad de viajes a la hora de caracterizar nodos clave del sistema de transporte urbano.

3 Construcción de la red multicapa

La idea es pasar de paradas individuales a zonas de intercambio (clusters de paradas cercanas), y luego construir una capa por modo de transporte.

3.1 1. Definición de zonas de intercambio (clustering espacial)

# Creamos objeto sf con coordenadas
transporte_sf <- transporte %>% 
  st_as_sf(coords = c("stop_lon", "stop_lat"), crs = 4326, remove = FALSE) %>% 
  st_transform(25830)  # sistema métrico (UTM zona 30)

coords_m <- transporte_sf %>% 
  st_coordinates() %>% 
  as.matrix()

# Clustering DBSCAN: paradas a menos de 150 m se agrupan
set.seed(123)
db <- dbscan(coords_m, eps = 150, minPts = 2)

transporte_sf$cluster_id <- db$cluster

# Tratamos el ruido (cluster 0) como clusters individuales
transporte_sf <- transporte_sf %>% 
  mutate(
    cluster_id = if_else(
      cluster_id == 0,
      max(cluster_id) + row_number(),
      cluster_id
    )
  )

# Zona de intercambio = centroide del cluster
zonas <- transporte_sf %>% 
  st_drop_geometry() %>% 
  group_by(cluster_id) %>% 
  summarise(
    stop_lat = mean(stop_lat),
    stop_lon = mean(stop_lon),
    # distrito asignado más frecuente en el cluster
    distrito_asignado = as.integer(names(which.max(table(distrito_asignado)))),
    modos_disponibles = paste(sort(unique(transport_mode)), collapse = ","),
    .groups = "drop"
  )

nrow(zonas)

## [1] 1558

head(zonas)

Cada cluster_id representa una zona de intercambio multimodal.

Representamos la distribución espacial aproximada (proyección geográfica simple):

zonas_sf <- zonas %>% 
  st_as_sf(coords = c("stop_lon", "stop_lat"), crs = 4326)

plot(zonas_sf["distrito_asignado"], main = "Zonas de intercambio por distrito")

3.2 2. Nodos por modo de transporte

modos <- sort(unique(transporte_sf$transport_mode))

zonas_modo <- transporte_sf %>% 
  st_drop_geometry() %>% 
  select(cluster_id, transport_mode) %>% 
  distinct() %>% 
  mutate(presente = 1) %>% 
  pivot_wider(
    names_from = transport_mode,
    values_from = presente,
    values_fill = 0
  )

zonas_attrs <- zonas %>% 
  left_join(zonas_modo, by = "cluster_id") %>% 
  mutate(across(all_of(modos), ~replace_na(., 0)))

kable(head(zonas_attrs), caption = "Atributos básicos de las zonas por modo")

Atributos básicos de las zonas por modo
cluster_id	stop_lat	stop_lon	distrito_asignado	modos_disponibles	bicimad	bus	metro
1	40.38368	-3.624114	2807918	bicimad,bus,metro	1	1	1
2	40.38096	-3.728060	2807911	bicimad,bus,metro	1	1	1
3	40.38626	-3.719795	2807912	bicimad,bus,metro	1	1	1
4	40.37886	-3.732537	2807911	bicimad,bus	1	1	0
5	40.40855	-3.672909	2807903	bicimad,bus	1	1	0
6	40.40380	-3.706547	2807901	bicimad,metro	1	0	1

3.3 3. Unimos la información de viajes por distrito

# Forzamos el mismo tipo en ambas columnas de distrito
zonas_attrs <- zonas_attrs %>% 
  mutate(distrito_asignado = as.integer(distrito_asignado))

viajes_distrito <- viajes_distrito %>% 
  mutate(distrito = as.integer(distrito))

# Join seguro
zonas_attrs <- zonas_attrs %>% 
  left_join(
    viajes_distrito,
    by = c("distrito_asignado" = "distrito")
  )

summary(zonas_attrs$viajes_medios_por_persona)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   1.956   2.004   2.038   2.041   2.063   2.135

En muchas zonas puede que no haya dato (si el distrito no aparece o hay problemas de asignación), pero para el análisis de red esto no es un problema; simplemente tendremos algunos NA.

3.4 4. Construcción de redes intra-capa

Para cada modo:

Consideramos solo las zonas donde ese modo está presente.
Unimos cada zona con sus k vecinos más cercanos (aquí k = 3) basándonos en la distancia euclídea de los centroides.
Asignamos un peso inversamente proporcional a la distancia.

# Centroides en coordenadas métricas
zonas_centroides <- zonas_sf %>% 
  st_transform(25830) %>% 
  st_coordinates() %>% 
  as_tibble() %>% 
  rename(x = X, y = Y) %>% 
  mutate(cluster_id = zonas$cluster_id)

build_edges_modo <- function(modo, k = 3) {
  zonas_modo_ids <- zonas_attrs %>% 
    filter(.data[[modo]] == 1) %>% 
    pull(cluster_id)

  if (length(zonas_modo_ids) <= 1) {
    return(tibble(modo = character(), from = integer(),
                  to = integer(), dist = numeric(), w = numeric()))
  }

  coords_modo <- zonas_centroides %>% 
    filter(cluster_id %in% zonas_modo_ids) %>% 
    arrange(cluster_id)

  mat <- as.matrix(coords_modo[, c("x", "y")])
  k_eff <- min(k, nrow(mat) - 1)
  knn <- FNN::get.knn(mat, k = k_eff)

  from <- rep(coords_modo$cluster_id, each = k_eff)
  to   <- coords_modo$cluster_id[c(knn$nn.index)]
  dist <- as.numeric(knn$nn.dist)

  tibble(
    modo  = modo,
    from  = from,
    to    = to,
    dist  = dist,
    w     = 1 / (1 + dist)
  ) %>% 
    distinct() %>% 
    filter(from != to)
}

edges_intra <- map_df(modos, build_edges_modo)

kable(head(edges_intra), caption = "Ejemplo de aristas intra-capa")

Ejemplo de aristas intra-capa
modo	from	to	dist	w
bicimad	1	288	649.2619	0.0015378
bicimad	1	300	321.5842	0.0031000
bicimad	1	12582	476.5779	0.0020939
bicimad	2	283	506.9860	0.0019686
bicimad	2	927	343.8273	0.0029000
bicimad	3	29	398.1358	0.0025054

3.5 5. Grafos por capa (lista `g_list`)

nodos_ids <- zonas_attrs$cluster_id

build_graph_modo <- function(modo) {
  e_modo <- edges_intra %>% 
    filter(modo == !!modo) %>% 
    select(from, to, w)

  g <- graph_from_data_frame(
    d = e_modo,
    directed = FALSE,
    vertices = zonas_attrs %>% 
      mutate(name = cluster_id)
  )

  missing_nodes <- setdiff(nodos_ids, as.integer(V(g)$name))
  if (length(missing_nodes) > 0) {
    g <- g + vertices(as.character(missing_nodes))
  }

  g
}

g_list <- map(modos, build_graph_modo)
names(g_list) <- modos

g_list

## $bicimad
## IGRAPH 6bce25f UN-- 1558 1401 -- 
## + attr: name (v/n), stop_lat (v/n), stop_lon (v/n), distrito_asignado
## | (v/n), modos_disponibles (v/c), bicimad (v/n), bus (v/n), metro
## | (v/n), parking (v/n), personas_totales (v/n), viajes_totales (v/n),
## | viajes_medios_por_persona (v/n), w (e/n)
## + edges from 6bce25f (vertex names):
##  [1]  1--  288  1--  300  1--12582  2--  283  2--  927  3--   29  3--  735
##  [8]  3--  248  4--  372  4-- 6182  4--15480  5--  583  5--  778  5--  626
## [15]  6--  296  6-- 1121  6--  813  7--  239  7--  327  7--  370  8--  209
## [22]  8--  152  8--  145 11-- 1182 11--  813 11-- 6031 12-- 6031 12--  115
## [29] 12-- 4280 13--   47 13--  174 13-- 9493 14-- 1256 14--  621 14--  814
## + ... omitted several edges
## 
## $bus
## IGRAPH 6bd1215 UN-- 1558 4036 -- 
## + attr: name (v/n), stop_lat (v/n), stop_lon (v/n), distrito_asignado
## | (v/n), modos_disponibles (v/c), bicimad (v/n), bus (v/n), metro
## | (v/n), parking (v/n), personas_totales (v/n), viajes_totales (v/n),
## | viajes_medios_por_persona (v/n), w (e/n)
## + edges from 6bd1215 (vertex names):
##  [1]  1--1289  1--1026  1--1303  2-- 652  2-- 283  2--  29  3-- 819  3-- 390
##  [9]  3-- 844  4-- 767  4--1288  4-- 638  5-- 158  5-- 422  5-- 969  7--  20
## [17]  7--  18  7--  18  8-- 997  8-- 785  8-- 464  9-- 734  9-- 365  9-- 588
## [25] 10--1121 10--1011 10-- 239 11--  31 11--1118 11--  33 12--  32 12-- 365
## [33] 12-- 532 13-- 296 13-- 938 13--1255 14-- 854 14-- 599 14--1182 15--  43
## + ... omitted several edges
## 
## $metro
## IGRAPH 6bd1ffc UN-- 1558 682 -- 
## + attr: name (v/n), stop_lat (v/n), stop_lon (v/n), distrito_asignado
## | (v/n), modos_disponibles (v/c), bicimad (v/n), bus (v/n), metro
## | (v/n), parking (v/n), personas_totales (v/n), viajes_totales (v/n),
## | viajes_medios_por_persona (v/n), w (e/n)
## + edges from 6bd1ffc (vertex names):
##  [1]  1--230  1--229  1--  2  2--239  2--239  2--248  3--372  3--112  3--187
## [10]  6--187  6--248  6--917  7--287  7-- 68  7--145 11--753 11--220 11-- 49
## [19] 12--235 12-- 47 12--174 13--249 13-- 82 13--423 17--814 17-- 59 17-- 57
## [28] 18-- 57 18--516 18--564 24-- 30 24--446 24--223 27--222 27--437 27--118
## [37] 28--348 28--523 28--243 30--170 30--145 30--573
## + ... omitted several edges
## 
## $parking
## IGRAPH 6bd2c18 UN-- 1558 161 -- 
## + attr: name (v/n), stop_lat (v/n), stop_lon (v/n), distrito_asignado
## | (v/n), modos_disponibles (v/c), bicimad (v/n), bus (v/n), metro
## | (v/n), parking (v/n), personas_totales (v/n), viajes_totales (v/n),
## | viajes_medios_por_persona (v/n), w (e/n)
## + edges from 6bd2c18 (vertex names):
##  [1]   7--   29   7--14152   7--  766   8--11008   8--  239   8--  245
##  [7]  27-- 1052  27--  814  27--  814  28--  446  28--  657  28--  164
## [13]  29--  573  29--  766  29--  941  41--  423  41--11008  41--  562
## [19]  29--   55  55--  423  55--  241  56--  939  56--  534  56--  678
## [25]  67--  613  67--14152  67--  630  71--  630  71--  814  71--  241
## + ... omitted several edges

4 Visualización multicapa con `muxViz`

En esta sección generamos la visualización clave: un gráfico 3D en el que se ven las capas (modos) apiladas, con los nodos alineados entre capas.

Para ello usamos la función plot_multiplex3D de muxViz. Se asume que muxViz está instalado; si no lo está, el código no fallará, pero mostrará un mensaje.

if (requireNamespace("muxViz", quietly = TRUE)) {
  library(muxViz)

  layer_colors <- c(
    bicimad = "#1a9641",
    bus     = "#d7191c",
    metro   = "#2c7bb6",
    parking = "#fdae61"
  )

  node_multimodalidad <- zonas_attrs %>% 
    mutate(
      modos_activos = rowSums(across(all_of(modos)))
    ) %>% 
    pull(modos_activos)

  node_size_values <- rescale(node_multimodalidad, to = c(0.5, 2))

  muxViz::plot_multiplex3D(
    g.list           = g_list,
    layer.colors     = layer_colors[names(g_list)],
    layer.labels     = "auto",     # <- clave: escalar, no vector
    node.size.values = node_size_values,
    node.size.scale  = 1.5,
    edge.size.scale  = 0.5,
    layout           = "fr",
    show.aggregate   = TRUE,
    show.nodeLabels  = FALSE
  )
} else {
  message("El paquete 'muxViz' no está instalado. Instálalo con devtools::install_github('manlius/muxViz') para ver la visualización 3D.")
}

Dado que la función plot_multiplex3D de muxViz presenta problemas en la instalación utilizada (bug interno con layer.labels), se ha implementado un gráfico multicapa alternativo usando tidygraph y ggraph, que representa las capas apiladas y las conexiones intra- e inter-capa de forma equivalente.

# Gráfico multicapa alternativo con tidygraph + ggraph

library(tidygraph)
library(ggraph)

# Orden de las capas y un índice para apilarlas verticalmente
layer_order <- modos
layer_index <- setNames(seq_along(layer_order), layer_order)

# Distancia vertical entre capas (en coordenadas proyectadas)
offset <- (max(zonas_centroides$y) - min(zonas_centroides$y)) * 1.2

# Nodos multilayer: mismo cluster_id repetido por capa, apilado en vertical
nodes_multi <- purrr::map_dfr(layer_order, function(m) {
  zonas_centroides %>%
    transmute(
      node       = paste(cluster_id, m, sep = "_"),  # id único por nodo-capa
      cluster_id = cluster_id,
      layer      = m,
      layer_idx  = layer_index[m],
      x          = x,
      y          = y + (layer_index[m] - 1) * offset
    )
})

# Aristas intra-capa: usamos edges_intra y las asociamos a cada capa
edges_intra_multi <- edges_intra %>%
  transmute(
    from  = paste(from, modo, sep = "_"),
    to    = paste(to,   modo, sep = "_"),
    type  = "intra",
    layer = modo
  )

# Aristas inter-capa: unen el MISMO cluster_id entre modos presentes (verticales)
inter_edges <- purrr::map_dfr(zonas_attrs$cluster_id, function(cid) {
  # modos presentes en esa zona
  fila <- zonas_attrs[zonas_attrs$cluster_id == cid, modos, drop = FALSE]
  modos_presentes <- modos[as.logical(unlist(fila))]
  
  if (length(modos_presentes) < 2) return(tibble())
  
  comb <- t(combn(modos_presentes, 2))
  tibble(
    from  = paste(cid, comb[, 1], sep = "_"),
    to    = paste(cid, comb[, 2], sep = "_"),
    type  = "inter",
    layer = "inter"
  )
})

edges_all <- bind_rows(edges_intra_multi, inter_edges)

# Grafo multilayer
g_multi <- tbl_graph(nodes = nodes_multi, edges = edges_all, directed = FALSE)

# Colores de capa (igual que antes)
layer_colors <- c(
  bicimad = "#1a9641",
  bus     = "#d7191c",
  metro   = "#2c7bb6",
  parking = "#fdae61"
)

ggraph(g_multi, layout = "manual", x = x, y = y) +
  geom_edge_link(aes(alpha = type, linetype = type), show.legend = TRUE) +
  geom_node_point(aes(color = layer), size = 1.2) +
  scale_edge_alpha_manual(values = c(intra = 0.3, inter = 0.8)) +
  scale_edge_linetype_manual(values = c(intra = "solid", inter = "dashed")) +
  scale_color_manual(values = layer_colors) +
  theme_void() +
  labs(
    title = "Red multicapa de transporte de Madrid",
    subtitle = "Capas (modos) apiladas verticalmente y conectadas por zonas de intercambio",
    color = "Capa",
    linetype = "Tipo de enlace",
    alpha = "Tipo de enlace"
  )

Este gráfico muestra claramente:

Las capas (bus, metro, bicimad, parking) como planos separados.
Los nodos (zonas de intercambio) alineados entre capas.
La estructura de conexiones en cada capa y la capa agregada.

# ===========================================
# Visualización 3D multicapa (inter-capa con scatter3d lines)
# ===========================================

if (requireNamespace("plotly", quietly = TRUE)) {

  library(plotly)
  library(dplyr)
  library(purrr)
  library(scales)
  library(tidyr)

  # ---- 1. Capas y eje Z
  layer_order <- modos
  layer_index <- setNames(seq_along(layer_order) - 1, layer_order)  # 0,1,2,3

  # ---- 2. Nodos multilayer: solo zonas donde ese modo existe
  nodes_multi_3d <- map_dfr(layer_order, function(m) {
    zonas_attrs %>%
      filter(.data[[m]] == 1) %>%                # solo zonas con ese modo
      select(cluster_id) %>%
      inner_join(zonas_centroides, by = "cluster_id") %>%
      transmute(
        node       = paste(cluster_id, m, sep = "_"),
        cluster_id = cluster_id,
        layer      = m,
        x_raw      = x,
        y_raw      = y
      )
  })

  # Normalizamos X, Y y asignamos Z
  nodes_multi_3d <- nodes_multi_3d %>%
    mutate(
      x = rescale(x_raw, to = c(0, 1)),
      y = rescale(y_raw, to = c(0, 1)),
      z = layer_index[layer]
    )

  # ---- 3. Multimodalidad por zona
  node_multimodalidad <- zonas_attrs %>%
    mutate(modos_activos = rowSums(across(all_of(modos)))) %>%
    select(cluster_id, modos_activos)

  # Clusters con al menos 2 modos
  clusters_interes <- node_multimodalidad %>%
    filter(modos_activos >= 2) %>%
    pull(cluster_id)

  # ---- 4. Aristas inter-capa: todos los pares de modos presentes
  edges_inter_3d <- map_dfr(clusters_interes, function(cid) {

    fila <- zonas_attrs[zonas_attrs$cluster_id == cid, modos, drop = FALSE]
    modos_presentes <- modos[as.logical(unlist(fila))]

    if (length(modos_presentes) < 2) return(tibble())

    comb <- t(combn(modos_presentes, 2))  # todos los pares posibles

    tibble(
      from       = paste(cid, comb[, 1], sep = "_"),
      to         = paste(cid, comb[, 2], sep = "_"),
      cluster_id = cid
    )
  })

  cat("Número de aristas inter-capa generadas:", nrow(edges_inter_3d), "\n")

  # ---- 5. Posiciones de nodos para las aristas
  pos_nodes <- nodes_multi_3d %>%
    select(node, x, y, z, layer)

  edges_plot <- edges_inter_3d %>%
    left_join(pos_nodes, by = c("from" = "node")) %>%
    rename(
      x  = x,
      y  = y,
      z  = z
    ) %>%
    left_join(pos_nodes, by = c("to" = "node"), suffix = c("", "_to")) %>%
    rename(
      xend = x_to,
      yend = y_to,
      zend = z_to
    ) %>%
    filter(
      !is.na(x), !is.na(y), !is.na(z),
      !is.na(xend), !is.na(yend), !is.na(zend)
    ) %>%
    select(x, y, z, xend, yend, zend, cluster_id)

  cat("Aristas con coordenadas válidas:", nrow(edges_plot), "\n")

  # ---- 6. Convertimos aristas a formato "long" para scatter3d lines
  # Para cada arista: (x,y,z) -> (xend,yend,zend) -> NA como separador
  if (nrow(edges_plot) > 0) {
    edges_long <- edges_plot %>%
      mutate(seg_id = row_number()) %>%
      mutate(
        x  = map2(x,  xend,  ~c(.x, .y, NA_real_)),
        y  = map2(y,  yend,  ~c(.x, .y, NA_real_)),
        z  = map2(z,  zend,  ~c(.x, .y, NA_real_))
      ) %>%
      select(seg_id, x, y, z) %>%
      unnest(c(x, y, z))
  } else {
    edges_long <- tibble(x = numeric(), y = numeric(), z = numeric())
  }

  cat("Puntos totales en edges_long:", nrow(edges_long), "\n")

  # ---- 7. Añadimos multimodalidad y tamaño de nodo
  nodes_multi_3d <- nodes_multi_3d %>%
    left_join(node_multimodalidad, by = "cluster_id") %>%
    mutate(
      modos_activos = replace_na(modos_activos, 1),
      size = rescale(modos_activos, to = c(4, 11))
    )

  # ---- 8. Colores por capa
  layer_colors <- c(
    bicimad = "#1a9641",
    bus     = "#d7191c",
    metro   = "#2c7bb6",
    parking = "#fdae61"
  )

  # ---- 9. Gráfico 3D
  p3d <- plot_ly()

  # Nodos
  p3d <- p3d %>%
    add_trace(
      data = nodes_multi_3d,
      type = "scatter3d",
      mode = "markers",
      x = ~x, y = ~y, z = ~z,
      color = ~layer,
      colors = layer_colors,
      marker = list(size = ~size),
      hoverinfo = "text",
      text = ~paste(
        "Zona:", cluster_id,
        "<br>Modo:", layer,
        "<br>Modos activos (total cluster):", modos_activos
      )
    )

  # Aristas inter-capa como líneas rojas gruesas
  if (nrow(edges_long) > 0) {
    p3d <- p3d %>%
      add_trace(
        data = edges_long,
        type = "scatter3d",
        mode = "lines",
        x = ~x, y = ~y, z = ~z,
        line = list(color = "black", width = 4),
        showlegend = FALSE,
        hoverinfo = "none"
      )
  } else {
    message("edges_long está vacío: no se han podido dibujar aristas inter-capa.")
  }

  p3d <- p3d %>%
    layout(
      title = "Red multicapa de transporte de Madrid (3D, enlaces entre capas)",
      scene = list(
        xaxis = list(title = "X (normalizada)"),
        yaxis = list(title = "Y (normalizada)"),
        zaxis = list(
          title = "Capa (modo de transporte)",
          tickvals = layer_index,
          ticktext = layer_order
        ),
        aspectmode = "cube"
      )
    )

  p3d

} else {
  message("Instala 'plotly' con install.packages('plotly') para ver la visualización 3D.")
}

## Número de aristas inter-capa generadas: 728 
## Aristas con coordenadas válidas: 728 
## Puntos totales en edges_long: 2184

Interpretación de la visualización 3D multicapa

La visualización tridimensional obtenida permite observar de forma explícita la estructura multicapa del sistema de transporte urbano de Madrid, donde cada plano horizontal representa un modo de transporte (bicimad, bus, metro y parking) y las conexiones verticales indican zonas de intercambio multimodal.

Separación clara de capas y estructura espacial

En primer lugar, se aprecia una separación nítida entre las cuatro capas, lo que confirma que cada modo de transporte presenta una distribución espacial propia:

La capa bicimad aparece más concentrada espacialmente, reflejando su carácter urbano y centralizado.
La capa bus muestra una cobertura más extensa, coherente con su función de red capilar que conecta barrios periféricos.
La capa metro se sitúa entre ambas, con una estructura más compacta que refleja la jerarquía de estaciones y líneas.
La capa parking aparece claramente diferenciada y menos densa, lo que concuerda con su función de infraestructura de apoyo más localizada.

Esta separación valida la decisión metodológica de modelar el sistema como una red multicapa, ya que una red agregada ocultaría estas diferencias estructurales.

Uniones verticales como indicador de multimodalidad

Las líneas verticales que conectan nodos entre capas representan zonas donde coexisten varios modos de transporte. Estas conexiones no aparecen de forma homogénea, sino que se concentran en un subconjunto reducido de zonas, lo que pone de manifiesto una fuerte heterogeneidad en la multimodalidad del sistema.

En términos de teoría de redes:

Estas zonas actúan como nodos altamente versátiles, ya que participan simultáneamente en varias capas.
Son puntos críticos para la integración funcional del sistema, permitiendo el transbordo eficiente entre modos.
Su posición estructural las convierte en puntos de control del flujo global de movilidad.

Visualmente, estas zonas destacan como “columnas” que atraviesan varias capas, reforzando la idea de que la conectividad global del sistema depende de un número limitado de nodos clave.

Versatilidad y roles multinodales

La figura ilustra de manera intuitiva el concepto de versatilidad de nodo:

La mayoría de los nodos aparecen en una sola capa, actuando como nodos especializados (por ejemplo, paradas de bus o estaciones de bici aisladas).
En contraste, un número reducido de nodos conecta tres o incluso cuatro capas, desempeñando un rol multinodal esencial.

Este patrón es consistente con la literatura sobre redes multicapas, donde se observa que la importancia de un nodo no puede evaluarse correctamente desde una sola capa. La visualización 3D confirma que los nodos versátiles no solo existen, sino que son estructuralmente visibles cuando se representan explícitamente las capas y sus interdependencias.

Implicaciones para vulnerabilidad y resiliencia

Desde el punto de vista de la vulnerabilidad del sistema, la visualización sugiere una conclusión relevante:

La conectividad multimodal del transporte madrileño depende de un conjunto relativamente pequeño de nodos altamente interconectados entre capas.
Una perturbación localizada en estos nodos (por ejemplo, cierre de un intercambiador o fallo simultáneo de varios modos en una zona) podría propagarse rápidamente al resto del sistema, afectando a múltiples capas.

Este comportamiento es coherente con los modelos de cascadas de fallo en redes interdependientes descritos en la literatura (Buldyrev et al., 2010; Duan et al., 2019), donde la interdependencia aumenta tanto la eficiencia como la fragilidad del sistema.

Valor añadido de la representación multicapa

Finalmente, esta visualización demuestra el valor añadido del enfoque multicapa frente a representaciones tradicionales:

Permite identificar de forma directa los intercambiadores críticos.
Hace visible la estructura de interdependencia entre modos.
Facilita la interpretación de conceptos abstractos como versatilidad, roles multinodales y vulnerabilidad sistémica.

En conjunto, el gráfico 3D no solo actúa como herramienta exploratoria, sino como una evidencia visual clara de que el sistema de transporte de Madrid funciona como una red interconectada de capas, cuya eficiencia y resiliencia dependen de la correcta gestión de sus nodos multimodales clave.

5 Centralidad y versatilidad de las zonas

5.1 1. Medidas de centralidad por capa

Calculamos grado y fuerza (suma de pesos) para cada zona en cada modo.

centralidad_por_capa <- map2_df(
  .x = g_list,
  .y = names(g_list),
  .f = function(g, modo) {
    tibble(
      cluster_id = as.integer(V(g)$name),
      modo       = modo,
      degree     = degree(g),
      strength   = strength(g, weights = E(g)$w)
    )
  }
)

kable(head(centralidad_por_capa), caption = "Centralidad (grado y fuerza) por capa")

Centralidad (grado y fuerza) por capa
cluster_id	modo	degree	strength
1	bicimad	5	0.0096077
2	bicimad	4	0.0090672
3	bicimad	4	0.0078289
4	bicimad	6	0.0134323
5	bicimad	7	0.0140216
6	bicimad	6	0.0139194

Distribución del grado por modo:

centralidad_por_capa %>% 
  ggplot(aes(x = degree, fill = modo)) +
  geom_histogram(alpha = 0.6, bins = 30, position = "identity") +
  facet_wrap(~ modo, scales = "free_y") +
  theme_minimal() +
  labs(
    title = "Distribución del grado por modo de transporte",
    x = "Grado",
    y = "Frecuencia"
  )

5.2 2. Versatilidad de grado

Definimos una medida de versatilidad basada en la entropía de Shannon de la distribución de grado entre modos.

mat_grado <- centralidad_por_capa %>% 
  select(cluster_id, modo, degree) %>% 
  pivot_wider(
    names_from = modo,
    values_from = degree,
    values_fill = 0
  ) %>% 
  arrange(cluster_id)

grado_mat <- as.matrix(mat_grado[, modos])
rownames(grado_mat) <- mat_grado$cluster_id

calc_versatilidad <- function(vec) {
  if (sum(vec) == 0) return(NA_real_)
  p <- vec / sum(vec)
  p <- p[p > 0]
  H <- -sum(p * log(p))
  H / log(length(vec))
}

versatilidad <- apply(grado_mat, 1, calc_versatilidad)

versatilidad_df <- tibble(
  cluster_id = as.integer(names(versatilidad)),
  versatilidad_grado = as.numeric(versatilidad)
) %>% 
  left_join(zonas_attrs, by = "cluster_id")

summary(versatilidad_df$versatilidad_grado)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0000  0.0000  0.0000  0.1477  0.0000  0.9983

5.2.1 Nodos más versátiles y más especializados

top_versatiles <- versatilidad_df %>% 
  arrange(desc(versatilidad_grado)) %>% 
  slice_head(n = 15) %>% 
  select(cluster_id, versatilidad_grado, modos_disponibles, distrito_asignado)

top_especializados <- versatilidad_df %>% 
  arrange(versatilidad_grado) %>% 
  slice_head(n = 15) %>% 
  select(cluster_id, versatilidad_grado, modos_disponibles, distrito_asignado)

kable(top_versatiles, caption = "Top 15 zonas más versátiles")

Top 15 zonas más versátiles
cluster_id	versatilidad_grado	modos_disponibles	distrito_asignado
164	0.9983119	bicimad,bus,metro,parking	2807907
446	0.9976183	bicimad,bus,metro,parking	2807904
101	0.9968869	bicimad,bus,metro,parking	2807907
562	0.9963891	bicimad,bus,metro,parking	2807920
41	0.9955380	bicimad,bus,metro,parking	2807907
56	0.9949672	bicimad,bus,metro,parking	2807904
283	0.9949672	bicimad,bus,metro,parking	2807903
239	0.9946468	bicimad,bus,metro,parking	2807901
814	0.9934804	bicimad,bus,metro,parking	2807904
71	0.9927376	bicimad,bus,metro,parking	2807904
339	0.9927376	bicimad,bus,metro,parking	2807911
534	0.9927376	bicimad,bus,metro,parking	2807908
7	0.9926388	bicimad,bus,metro,parking	2807901
251	0.9892423	bicimad,bus,metro,parking	2807916
116	0.9863237	bicimad,bus,metro,parking	2807907

kable(top_especializados, caption = "Top 15 zonas más especializadas")

Top 15 zonas más especializadas
cluster_id	modos_disponibles	distrito_asignado
9	bus	2807903
10	bus	2807909
15	bus	2807921
17	metro	2807921
19	bus	2807921
20	bus	2807921
21	bus	2807921
22	bus	2807908
26	bus	2807911
31	bus	2807921
32	bus	2807921
33	bus	2807921
34	bus	2807911
36	bus	2807911
38	bus	2807920

Interpretación cualitativa (a redactar en la memoria):

Zonas con alta versatilidad: concentran paradas de varios modos y son relevantes en varias capas.
Zonas especializadas: la conectividad se concentra en un único modo.

5.3 3. Relación entre versatilidad e intensidad de viajes

versatilidad_df %>% 
  ggplot(aes(x = viajes_medios_por_persona, y = versatilidad_grado)) +
  geom_point(alpha = 0.4) +
  geom_smooth(method = "loess", se = FALSE) +
  theme_minimal() +
  labs(
    title = "Versatilidad de las zonas vs viajes medios por persona",
    x = "Viajes medios por persona (distrito)",
    y = "Versatilidad de grado"
  )

Comprobamos cuántas observaciones completas hay y calculamos la correlación solo si tiene sentido estadísticamente (al menos 3 observaciones completas):

datos_cor <- versatilidad_df %>% 
  filter(!is.na(versatilidad_grado),
         !is.na(viajes_medios_por_persona))

n_obs <- nrow(datos_cor)
n_obs

## [1] 1558

if (n_obs >= 3) {
  cor_test <- cor.test(
    datos_cor$versatilidad_grado,
    datos_cor$viajes_medios_por_persona
  )
  cor_test
} else {
  cat("No hay suficientes observaciones completas para estimar de forma fiable la correlación entre versatilidad y viajes.
")
}

## 
##  Pearson's product-moment correlation
## 
## data:  datos_cor$versatilidad_grado and datos_cor$viajes_medios_por_persona
## t = 0.13627, df = 1556, p-value = 0.8916
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.04621554  0.05310757
## sample estimates:
##         cor 
## 0.003454533

Esto evita el error de “not enough finite observations” y hace el análisis más robusto.

6 Vulnerabilidad de la red de transporte

En lugar de implementar un modelo muy complejo de cascadas, realizamos un análisis sencillo pero informativo de robustez estructural:

Tomamos la red agregada (unión de todas las capas).
Eliminamos nodos de forma dirigida (los más centrales) en una capa y observamos cómo disminuye el tamaño de la componente gigante agregada.

6.1 1. Red agregada

# Matrices de adyacencia por capa
adj_list <- lapply(g_list, function(g) {
  as.matrix(as_adj(g, attr = "w", sparse = FALSE))
})

# Aseguramos mismo orden de nodos
orden_nodos <- order(as.integer(V(g_list[[1]])$name))
adj_list <- lapply(adj_list, function(m) m[orden_nodos, orden_nodos])

adj_agregada <- Reduce("+", adj_list)

g_agregado <- graph_from_adjacency_matrix(adj_agregada, mode = "undirected", weighted = TRUE)

comp <- components(g_agregado)
tamano_gcc_inicial <- max(comp$csize)
tamano_gcc_inicial

## [1] 1558

6.2 2. Robustez frente a fallos en una capa

Analizamos qué ocurre si vamos eliminando nodos con mayor grado en la capa metro y miramos el tamaño relativo de la componente gigante en la red agregada.

g_metro <- g_list[["metro"]]

deg_metro <- degree(g_metro)
orden_fallo <- names(sort(deg_metro, decreasing = TRUE))  # orden de fallo: más conectados primero

fracciones <- seq(0, 0.5, by = 0.05)

resultado_robustez <- map_df(fracciones, function(f) {
  n_remove <- ceiling(f * length(orden_fallo))
  nodos_fuera <- orden_fallo[seq_len(n_remove)]

  g_agregado_reducido <- delete_vertices(g_agregado, nodos_fuera)
  comp_red <- components(g_agregado_reducido)
  gcc_red <- if (length(comp_red$csize) == 0) 0 else max(comp_red$csize)

  tibble(
    frac_eliminada = f,
    tamano_gcc_rel = gcc_red / tamano_gcc_inicial
  )
})

resultado_robustez %>% 
  ggplot(aes(x = frac_eliminada, y = tamano_gcc_rel)) +
  geom_line() +
  geom_point() +
  theme_minimal() +
  scale_y_continuous(labels = percent_format(accuracy = 1)) +
  scale_x_continuous(labels = percent_format(accuracy = 1)) +
  labs(
    title = "Robustez de la red agregada ante fallos en la capa METRO",
    x = "Fracción de nodos eliminados en METRO",
    y = "Tamaño relativo de la componente gigante agregada"
  )

Una caída rápida del tamaño de la componente gigante indica alta vulnerabilidad ante fallos dirigidos en nodos centrales de la capa metro.

7 Conclusiones

En esta memoria se ha construido y analizado una red multicapas de transporte para Madrid utilizando exclusivamente los datos proporcionados:

Las paradas se han agrupado en zonas de intercambio mediante clustering espacial.
Se ha definido una red multicapa con capas para bus, metro, bicimad y parking.
Se ha calculado la centralidad de las zonas en cada capa y un indicador de versatilidad de grado que muestra qué zonas son importantes en varios modos a la vez.
Se ha generado una visualización multicapa con muxViz que hace explícita la estructura en capas.
Se ha estudiado la robustez estructural de la red agregada frente a fallos dirigidos en nodos centrales de la capa metro.

A partir de aquí, se podrían ampliar los análisis:

Incorporando información de líneas reales (no solo proximidad espacial).
Usando otras medidas de centralidad (betweenness, eigenvector, PageRank).
Implementando modelos de cascadas más cercanos a la literatura de redes interdependientes.

El código está preparado para ejecutarse directamente (asumiendo que los CSV y los paquetes necesarios están disponibles) y centrado únicamente en los datos reales de Madrid, tal y como se pedía en el enunciado.

Redes Multicapas en el Transporte de Madrid

Tu nombre

25/01/2026

1 Introducción

2 Carga y exploración de datos

2.1 Lectura de los ficheros

2.2 Descripción básica de `transporte_madrid_consolidado`

2.3 Descripción básica de `viajes_madrid`

2.3.1 Indicador de intensidad de viajes por distrito

2.3.2 Indicador de viajes medios diarios por persona y distrito

3 Construcción de la red multicapa

3.1 1. Definición de zonas de intercambio (clustering espacial)

3.2 2. Nodos por modo de transporte

3.3 3. Unimos la información de viajes por distrito

3.4 4. Construcción de redes intra-capa

3.5 5. Grafos por capa (lista `g_list`)

4 Visualización multicapa con `muxViz`

5 Centralidad y versatilidad de las zonas

5.1 1. Medidas de centralidad por capa

5.2 2. Versatilidad de grado

5.2.1 Nodos más versátiles y más especializados

5.3 3. Relación entre versatilidad e intensidad de viajes

6 Vulnerabilidad de la red de transporte

6.1 1. Red agregada

6.2 2. Robustez frente a fallos en una capa

7 Conclusiones

Redes Multicapas en el Transporte de Madrid

Tu nombre

25/01/2026

1 Introducción

2 Carga y exploración de datos

2.1 Lectura de los ficheros

2.2 Descripción básica de transporte_madrid_consolidado

2.3 Descripción básica de viajes_madrid

2.3.1 Indicador de intensidad de viajes por distrito

2.3.2 Indicador de viajes medios diarios por persona y distrito

3 Construcción de la red multicapa

3.1 1. Definición de zonas de intercambio (clustering espacial)

3.2 2. Nodos por modo de transporte

3.3 3. Unimos la información de viajes por distrito

3.4 4. Construcción de redes intra-capa

3.5 5. Grafos por capa (lista g_list)

4 Visualización multicapa con muxViz

5 Centralidad y versatilidad de las zonas

5.1 1. Medidas de centralidad por capa

5.2 2. Versatilidad de grado

5.2.1 Nodos más versátiles y más especializados

5.3 3. Relación entre versatilidad e intensidad de viajes

6 Vulnerabilidad de la red de transporte

6.1 1. Red agregada

6.2 2. Robustez frente a fallos en una capa

7 Conclusiones

2.2 Descripción básica de `transporte_madrid_consolidado`

2.3 Descripción básica de `viajes_madrid`

3.5 5. Grafos por capa (lista `g_list`)

4 Visualización multicapa con `muxViz`